
function dx = altitude_plant(t, x)
    % 状态解码
    phi = x(1:3); % 姿态角
    phi_dot = x(4:6); % 姿态角速度
    xf = x(7:9);
    xf_dot = x(10:12);
    gamma_hat = x(13); % 自适应参数
    k_hat = x(14);
    %% 仿真参数
    % 系统参数 J C
    I11 = 1.25; I22 = 1.25; I33 = 2.5; % kg·m^2
    [J, C] = computeQuadrotorMatrices(phi(1), phi(2), phi(3), phi_dot(1), ...
                                        phi_dot(2), phi_dot(3), I11, I22, I33);
    
    % 控制参数
    T_sigma1 = 1; T_sigma2 = 1; T_s = 4; T_lambda = 4;
    sigma1 = 5/3; sigma2 = 3/2; sigma3 = 37/39; sigma4 = 39/37;
    lambda1 = 1/2; lambda2 = 3/2;
    pi_0 = 0.1; eta_d = 0.9; eta_lambda = 0.9;
    cd = 84; c_lambda = 3;  g = 3/2; 
    
    alpha_s2 = 1/( eta_lambda * T_lambda * (1-lambda1) * 2^(lambda1/2-3/2) );
    beta_s2 = (2^(lambda2/2+3/2)) / ( eta_lambda * T_lambda * (lambda2 - 1) );
    
    alpha_d2 = 1/( eta_d * T_sigma2 * (1-sigma3) * 2^(sigma3/2-3/2) );
    beta_d2 = (2^(sigma4/2+3/2)) / ( eta_d * T_sigma2 * (sigma4 - 1) );

    % 期望姿态和外部扰动
    phi_d       = @(t) 0.1 * [sin(t); cos(t); 0]; % 期望姿态 (rad)
    phi_d_dot   = @(t) 0.1 * [cos(t); -sin(t); 0]; % 期望角速度
    phi_d_ddot  = @(t) -0.1 * [sin(t); cos(t); 0]; % 期望角加速度
    d       = @(t) 0.1 * [sin(t); cos(t); sin(2*t)]; % 外部扰动 (N·m)

    % 期望值
    phi_d_k = phi_d(t);
    phi_d_dot_k = phi_d_dot(t);
    phi_d_ddot_k = phi_d_ddot(t);
    tao_d = J\d(t);

    % 误差计算
    phi_e = phi - phi_d_k;
    phi_e_dot = phi_dot - phi_d_dot_k;
    
    % 辅助变量
    e1 = phi_e - xf;
    e2 = phi_e_dot - xf_dot;
    % 滑模变量
    Ke = diag( [k(e1(1)), k(e1(2)), k(e1(3))] ) ;
    Ke_tail = diag( [k_tail(e1(1)), k_tail(e1(2)), k_tail(e1(3))] ) ;

    Kx1 = diag( [k(phi_e(1)), k(phi_e(2)), k(phi_e(3))] ) ;
    Kx1_tail = diag( [k_tail(phi_e(1)), k_tail(phi_e(2)), k_tail(phi_e(3))] ) ;

    s = Ke * e1 + sig(e2,sigma2);

    % DO
    tao_d_hat = (1/sigma2)*(Ke_tail + Ke) * sig(e2,2-sigma2) + k_hat * sign(s) ...
    + (1/sigma2)*diag(( norm(e2)^(1-sigma2) ) * ones(1, 3)) * ( alpha_d2 * sig(s,sigma3) + beta_d2 * sig(s,sigma4) );
    
    % 控制律
    tao = C * phi_dot + J * phi_d_ddot_k - J * tao_d_hat - gamma_hat * J * sign(s) ...
    - 1/g*( Kx1+Kx1_tail )*sig(phi_e_dot,2-g) - 1/g * J * diag(( norm(phi_e_dot)^(1-g) ) * ones(1, 3)) * ( alpha_s2 * sig(s,lambda1) + beta_s2 * sig(s,lambda2) );

    % 干扰估计更新
    % tau_hat_dot = cd * (abs(s) - eta_d * tau_hat);

    % 自适应参数更新
    gamma_hat_dot = g * norm(s) * norm(phi_dot)^(g-1) - c_lambda * gamma_hat;

    k_hat_dot = sigma2 * norm(s) * norm(e2)^(sigma2 - 1) - cd * k_hat;

    % 动力学更新
    phi_ddot = -J \ C * phi_dot + J \ tao + tao_d - phi_d_ddot_k;

    %辅助变量更新
    xf_ddot = -J \ C * phi_dot + J \ tao + tao_d_hat - phi_d_ddot_k;

    % 状态导数
    dx = [phi_dot; phi_ddot; xf_dot; xf_ddot; gamma_hat_dot; k_hat_dot];
end
